home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / util.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-08-22  |  3.6 KB  |  138 lines
  1. /* integration/util.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. static inline
  21. void update (gsl_integration_workspace * workspace,
  22.          double a1, double b1, double area1, double error1,
  23.          double a2, double b2, double area2, double error2);
  24.  
  25. static inline void
  26. retrieve (const gsl_integration_workspace * workspace, 
  27.       double * a, double * b, double * r, double * e);
  28.  
  29.  
  30.  
  31. static inline
  32. void update (gsl_integration_workspace * workspace,
  33.          double a1, double b1, double area1, double error1,
  34.          double a2, double b2, double area2, double error2)
  35. {
  36.   double * alist = workspace->alist ;
  37.   double * blist = workspace->blist ;
  38.   double * rlist = workspace->rlist ;
  39.   double * elist = workspace->elist ;
  40.   size_t * level = workspace->level ;
  41.  
  42.   const size_t i_max = workspace->i ;
  43.   const size_t i_new = workspace->size ;
  44.  
  45.   const size_t new_level = workspace->level[i_max] + 1;
  46.  
  47.   /* append the newly-created intervals to the list */
  48.   
  49.   if (error2 > error1)
  50.     {
  51.       alist[i_max] = a2;    /* blist[maxerr] is already == b2 */
  52.       rlist[i_max] = area2;
  53.       elist[i_max] = error2;
  54.       level[i_max] = new_level;
  55.       
  56.       alist[i_new] = a1;
  57.       blist[i_new] = b1;
  58.       rlist[i_new] = area1;
  59.       elist[i_new] = error1;
  60.       level[i_new] = new_level;
  61.     }
  62.   else
  63.     {
  64.       blist[i_max] = b1;    /* alist[maxerr] is already == a1 */
  65.       rlist[i_max] = area1;
  66.       elist[i_max] = error1;
  67.       level[i_max] = new_level;
  68.       
  69.       alist[i_new] = a2;
  70.       blist[i_new] = b2;
  71.       rlist[i_new] = area2;
  72.       elist[i_new] = error2;
  73.       level[i_new] = new_level;
  74.     }
  75.   
  76.   workspace->size++;
  77.  
  78.   if (new_level > workspace->maximum_level)
  79.     {
  80.       workspace->maximum_level = new_level;
  81.     }
  82.  
  83.   qpsrt (workspace) ;
  84. }
  85.  
  86. static inline void
  87. retrieve (const gsl_integration_workspace * workspace, 
  88.       double * a, double * b, double * r, double * e)
  89. {
  90.   const size_t i = workspace->i;
  91.   double * alist = workspace->alist;
  92.   double * blist = workspace->blist;
  93.   double * rlist = workspace->rlist;
  94.   double * elist = workspace->elist;
  95.  
  96.   *a = alist[i] ;
  97.   *b = blist[i] ;
  98.   *r = rlist[i] ;
  99.   *e = elist[i] ;
  100. }
  101.  
  102. static inline double
  103. sum_results (const gsl_integration_workspace * workspace);
  104.  
  105. static inline double
  106. sum_results (const gsl_integration_workspace * workspace)
  107. {
  108.   const double * const rlist = workspace->rlist ;
  109.   const size_t n = workspace->size;
  110.  
  111.   size_t k;
  112.   double result_sum = 0;
  113.  
  114.   for (k = 0; k < n; k++)
  115.     {
  116.       result_sum += rlist[k];
  117.     }
  118.   
  119.   return result_sum;
  120. }
  121.  
  122. static inline int
  123. subinterval_too_small (double a1, double a2, double b2);
  124.  
  125. static inline int
  126. subinterval_too_small (double a1, double a2, double b2)
  127. {
  128.   const double e = GSL_DBL_EPSILON;
  129.   const double u = GSL_DBL_MIN;
  130.  
  131.   double tmp = (1 + 100 * e) * (fabs (a2) + 1000 * u);
  132.  
  133.   int status = fabs (a1) <= tmp && fabs (b2) <= tmp;
  134.  
  135.   return status;
  136. }
  137.  
  138.